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ABSTRACT 



A three-dimensional (3D) coupled normal mode model for studying 
sound propagation in a complex coastal environment is developed. 
This development corresponds to a significant upgrade of an earlier 
version of the model in which a flat, rigid bottom was used. By 
imposing the general boundary conditions for an irregular, non- 
rigid bottom, the coupling coefficient integrals in the system of 
differential equations governing the mode amplitude are re- 
formulated. The model upgrade entails a numerical implementation of 
the revised formulae. With the improved physics, this latest 
version is capable of modeling the 3D acoustic wave-field in 
shallow water where sound speed, water depth and sediment 
properties can vary with horizontal location. To demonstrate this 
enhanced capability, the model is used here to simulate the 
interactions of the normal modes as they propagate up a sloping 
bottom. 
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I . INTRODUCTION 

A. BACKGROUND 

There are three approaches to model three-dimensional (3D) 
sound propagation in the ocean: ray theory, parabolic 
equation approximation and normal mode theory. 

Ray theory gives an approximate, planewave- like solution 
to the wave equation, which is valid at high enough 
frequencies and in media with gradual variations. The ray 
solution is constructed by raytracing. The acoustic rays 
provide for a visual, physical description of sound 
transmission in the ocean. The ray solution, however, 
neglects sound diffraction and dispersion and thus needs 
corrections near caustics and turning points. These 
corrections may sometimes be mathematically complicated. The 
Hamiltonian Acoustic Ray Tracing Program for the Ocean (HARPO) 
is the only 3D ray theory model available today. This computer 
code was originally developed by Jones et al . [Ref. 1] for the 
computation of 3D rays. 

The parabolic equation approximation method (PE) was 
introduced by Tappert [Ref. 2] . PE is a "full -wave" method 
that accounts for both sound diffraction and dispersion. It 
provides for numerical solutions to the wave equation which 
are accurate for energy propagating at low grazing angles. The 
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accuracy generally degrades as the angle increases. The 
backscattered energy is generally neglected in this 
approximation. A versatile 3D PE model has been developed by 
Lee et al . [Ref. 3] using an implicit finite difference 
scheme. Another 3D PE model was developed earlier by Baer 
[Ref. 4] which uses a split-step Fourier algorithm. The PE 
model of Lee has a wider angle capability, i.e., it models 
sound energy travelling at steeper angles more accurately. 

Finally, normal mode theory describes sound propagation as 
a collection of eigenfunctions, called normal modes, which are 
a natural set of vertical vibration modes in the sound 
channel. Just like PE, normal mode theory is a "full -wave" 
approach. The original normal mode theory assumes a 
horizontally stratified propagation medium. This assumption is 
valid for many short-range, deep-water cases, where range and 
azimuthal variations are negligibly small. Pierce [Ref. 5] 
extended the theory to account for horizontal sound speed, 
bathymetry and bottom-property variations. These variations 
produce mode coupling phenomena (in which energy exchange 
between modes takes place) . A 3D coupled normal mode model has 
been developed by Chiu and Ehret [Ref. 6] . This model is 
capable of simulating mode-mode interactions due to a 3D 
varying sound speed field. The effects of bottom bathymetry 
variations and sediment property, however, are not modeled. 
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B. THESIS OBJECTIVES AND OUTLINE 



The main objective of this thesis is to improve the Chiu- 
Ehret [Ref. 6] 3D coupled normal mode model by including the 
effects of bathymetry variations and sediment properties on 
sound propagation. The improved model is useful for studying 
sound propagation in shallow water environments where 
significant bottom interaction is expected. 

In Chapter II, 3D coupled mode theory is first reviewed. 
The formulae for the mode coupling coefficients in the system 
of differential equations governing the mode amplitude 
functions are derived. In the derivation, the general boundary 
conditions for an irregular, non-rigid bottom are used. 

In Chapter III, alternative expressions for the mode 
coupling coefficients are derived. These expressions allow for 
an easier numerical implementation. The improved model is used 
to examine the effects of a sloping bottom on upslope sound 
propagation. The validity of the adiabatic approximation is 
also examined. Conclusions are given in Chapter IV. 

The products coming out of this thesis are computer 
subroutines to include bathymetry variations and sediment 
properties in the 3D coupled mode model of Chiu and Ehret 
[Ref. 6] . The new routines are listed in the Appendix. 
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II. 



3D COUPLED NORMAL MODE THEORY 



In the mathematical formulation that follows, a 
cylindrical coordinate system will be used (see Fig. 1) . The 
z-axis is perpendicular to the ocean surface and is positive 
downward, r is range from the source location (i.e., the 
origin) and 8 is the azimuthal angle (positive clockwise) . 
Sound speed in the water column, c, , is a function of r, z and 
8, where sound speed in the sediment, c 2 , is assumed to be r 
and 8 dependent only. The density of the water column, p 1# is 
considered to be constant. The density of the sediment, p 2 , is 
also considered to be constant. The water- sediment interface 
is located at z=H(r,0). 

A. THE MATHEMATICAL PROBLEM 

In the case of isodensity layers, the 3D, homogeneous, 
monofrequency Helmholtz Equation governing the acoustic 
pressure, p, is: 

^piz , r ,0) + k 2 (z, i, 0)p(z, r, 0) = 0 (D 

where k (z, r, 8) =co/c (z, r, 8) is the acoustic wavenumber, w is the 
source angular frequency and c is sound speed (C[ in the water 
layer and c 2 in the sediment layer) . 

A quasi - separable solution to Eq. (1) is postulated, which 
is locally a linear combination of normal modes or depth 
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Figure 1. The model geometry 
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functions, Z n : 



p(z,r,0) = ]T R n (r, 0) Z n (z; r, 0) 



(2) 



n 



where R n are the mode amplitude functions and n is the mode 
number. 

The normal modes Z n are required to satisfy the depth 
equation: 



where k n is the horizontal component of the wavenumber 
(eigenvalue) associated with the n* mode. 

It can be easily shown, using the boundary conditions for 
Z n (to be derived next) and the depth equation (Eq. (3) ) , that 
the normal modes form a complete set of orthogonal functions, 
with the inverse of the medium density p as a weighting 
function in the normalization: 



where 6 nm is the Kronecker delta. Note that the integration is 
carried over the entire depth from 0 to <». Also, note that the 



dz, 




(3) 



OO 




(4) 
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density p in the model is considered to be constant in each 
layer (p, in the water and p 2 in the sediment) . 

B. THE BOUNDARY CONDITIONS 

The appropriate boundary condition for the acoustic 
pressure at the sea surface is 

p 1 ( z=0 ; r , 0) =0 (5) 

where the subscript 1 denotes the water column. This pressure 
release condition implies that the normal modes Z n must also 
be zero at the sea surface, i.e., 

Z n ( z-0 ; r , 6) = 0 (6) 

It also implies that the horizontal derivatives of Z n at z=0 
are zero, i.e., 

dZjz= Or, 9) = Q (?) 

dr 

1 dZjz= 0;r,6) = (8) 

r 80 

At the interface between the sediment and the water 
column, i.e., at z=H(r,0), the boundary conditions are 
continuity of pressure and continuity of the particle velocity 
component normal to the interface: 
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p 1 (z=H;r,Q) = p 2 ( z=H; r , 0) 



(9) 



1 ' A 1 A 
— Vp 1 . n - — Vp. . n 

Pi P 2 



( 10 ) 



The subscripts 1 and 2 denote the water column and the 
sediment, respectively. 

The unit directional vector n normal to the bottom 
interface is 



A dHU, 9) A 1 3tf(r,9) { 

£ = V(z-tf(r,6) ) = dr r 30 

IV ( z-tf (r , 0) I [ 1 |( 3//(r,0) l dHU, 6) )2 p 

3r r 30 



( 11 ) 



where z,r and 6 are the unit directional vectors associated 
with the z, r and 6 directions, respectively. In the case of 
a small bottom slope, the boundary condition of Eq. (10) can 
be approximated by 

1 3p t (z=tf; r, 9) _ i 3p 2 (z=tf;r,0) ( 
p x dz p 2 3z 



The small slope approximation is accurate when 
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(13) 



\ dHU^l \ < 1 , |1 dH ^' d) l < 1 

dr r d0 



Following Eq. (9) and Eq. (12), we obtain the following 
boundary conditions at the water- sediment interface for the 
normal modes : 



1 dZ ln (z=H;r,Q) = i dZ 2n (z=/f; r , 9) 
p x dz p 2 dz 



Z ln ( z=H; r , 0) = Z 2n ( z=tf, r , 0) (15) 

These boundary conditions hold for each individual normal mode 
because they are orthogonal functions. 

The boundary condition for p at z -* oo is 

p 2 ( z-<», r, 0) = 0 (16) 



This implies that the ’normal modes and their horizontal 
derivatives are also zero as z -* oo : 

Z 2n (z~°o;r,0) =0 (17) 



dZ 2n (z~°°; r,6) 
dr 



= 0 



(18) 
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(19) 



1 dZ 2n {z~°°; r ,0) 

r ae 



= o 



C. MODE COUPLING COEFFICIENTS DERIVATION 

Substituting Eq. (2) into Eq. (1), multiplying by 
Z m (z;r ,9)/p, integrating over the entire depth and finally 
rearranging terms, we obtain the coupled mode equations 
governing the mode amplitude functions: 



where the two mode coupling coefficients are defined as 




n 



( 20 ) 



00 




0 



( 21 ) 



and 



oo 




0 



( 22 ) 



V h is the horizontal gradient operator, i.e.. 
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(23) 



A g 



A 

+ 0 - 



_d_ 

00 



The mode coupling coefficients are measures of the 
significance of exchange of acoustic energy between modes 
resulting from horizontal variations in the medium. As the 
variations become stronger, the coupling coefficients become 
larger and so is the energy exchange. In the case of a 
completely range independent medium, the coupling 
coefficients are identically zero and the RHS of Eq. (20) 
vanishes. In such case, the modes propagate independently of 
each other. For range -dependent cases, the neglect of mode 
coupling leads to the adiabatic approximation [Ref. 5] . 

Cylindrical spreading can be removed from the coupled mode 
equation (Eq. (20) ) by replacing the mode amplitude function 
R n (r,0) with P n (r , 6) /r I/2 . The result is 



[^ + V(r,e) + 4^ + -U]P„(r, 0) - 

dr 2 r 2 ad 2 4 r 2 

= £ ^(r,©) . [V A P n (r , 0) -r Pn ^' Q) ] +B an (r,e) Pjr.B)) 



(24) 
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III. THE NUMERICAL MODEL AND EXAMPLE RUNS 



In this chapter, the procedure to upgrade the Chiu-Ehret 
model [Ref. 6] is discussed. The upgrade has entailed the 
derivation of alternative expressions for the mode coupling 
coefficients and the generation of new code to compute these 
coefficients based on the alternative expressions. 

The numerical results from two simple example model runs 
associated with two different bottom slopes are also presented 
in this chapter. Both cases deal with upslope propagation in 
isospeed layers. These runs have allowed for an examination of 
the coupling between modes caused by bathymetry change. In 
addition, they have allowed for an examination of the validity 
of the adiabatic approximation. 



A. THE CHIU-EHRET APPROACH 

In the far field, i.e., kr>>l , the coupled mode equation 
(Eq. (24)) for the mode amplitude functions, can be recast as 



[ — +ic 2 (r,0) + — — ] P n (r,0) i 
dr 2 m r 2 50 2 n 



EG„(r,6) .V t P„(r,e).B m (r,6)P„(r,0)) 



( 25 ) 



In the Chiu-Ehret model [Ref. 6], P n is decomposed as 
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P n (r, 0) = U n {r, 0) e 



j<t>„(r,e) 



r 

4> n ( r , 0 ) = J k n (r,Q) dr 
0 



( 26 ) 



where U n is the slowly varying complex envelope of P n 
modulating the rapidly varying two-dimensional (2D) adiabatic 
solution, i.e., exp(j$ n ), and <p n is the adiabatic phase. The 
Chiu-Ehret model numerically computes the envelopes U n instead 
of P n using Runge-Kutta schemes. 

B. ALTERNATIVE EXPRESSIONS FOR THE MODE COUPLING COEFFICIENTS 

For simpler numerical implementation, the expressions of 
the mode coupling coefficients in Eq. (21) and Eq. (22) are 
rewritten in alternative forms. These alternative forms do not 
require integrations of expressions involving the horizontal 
derivatives of normal modes. In the following, the derivation 
of these alternative forms is presented. 

1. VECTOR MODE COUPLING COEFFICIENT, ^ 
a. Case of m£n 

Applying the horizontal gradient operator V h to 
both sides of the depth equation (Eq. (3)), we get 
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+ 2 [k(z,r,Q) V h k(z,r,Q) - k n (r, 0) V h k n (r , 0) ] Z n 



( 27 ) 



dz 2 



[k 2 ( z, r , 0) - ic 2 n (r , 0) ] = 0 



Multiplying Eq. (27) by Z m (z;r ,6)/p and then integrating over 
the entire depth, we get 



oo r32 n 7 00 

T Az £-22 dz + f (k 2 ~k 2 ) —Z n V h Z n dz = 

J p m dz 2 J n p n h n 



2K - / | 

o p 



^ dz 



( 28 ) 



In order to recast the first term of Eq. (28) in a form useful 
for this derivation, we first use integration by parts twice 
with respect to z. The resulting expression, after some 
lengthy manipulations, is 



f-i 

J D 



1 „ &\Z n 



Z- 2_£ dz = f — 

p m dz 2 ip 



V,z n dz * 

h n dz 2 Pi lm 32 0 



\H(r,Q) 



— Z. 

P 2 



dw h z. 



2m 



h^2n 

dz 



00 


\H(r, 0) 


_ i dz lm 

_ w h^ln 


1 dz 

v h ^ 2n 



H(r,Q) 



( 29 ) 
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Again, the subscripts 1 and 2 denote the water column and the 
sediment, respectively. 

Application of the boundary conditions Eqs . (6), (7) , (8), 
(14), (15), (17), (18) and (19) to Eq. (29), yields 
subsequently 



Replacing the first term of Eq. (28) by Eq. (30) and then 
making use of the depth equation (Eq. (3)), we obtain the 
following alternative expression for the vector mode coupling 
coefficient, for m£n : 





( 30 ) 





z=H(z,Q) 
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2 



A A 

Ann( r ' 0 > = P wn r + Ymn 0 = 



1 . ,dV h Z ln dV h Z 2ns 

Z 1 m' ^ / 



dz 



dz 



k 2 {z,d)-k 2 { 

dZ 



[2 f— 

r , 0) J P 



* V /A A, z n dz + 



z=//(r, 0) 



“ <— V u -tW 



0Z p. 



_1 

P 2 



Z=H{Z , 6 ) ^ 

( 31 ) 



or equivalently, in light of the boundary conditions Eq. (14) 
and Eq. (15) , 



Ann ( r ' 9 ) 

1_ 

Pi 





A 

Y inn 0 = 


2 


oo 

- [2 {l kV h k An Z n dz + 
o M 


k n 2 ( r , 0) -k m 2 (r, 0) 


. 0 Vz lin 
lm dz 


z.»',,e, < pi p2 > 


^ Z lm y r 1 

g z v h^ln z-H(r,6) J 



( 32 ) 

The above expression only involves Z n and not their horizontal 
derivatives in the integrands. Therefore, the corresponding 
numerical evaluations are more efficient. 

The last two terms of Eq. (32) express the direct 
contribution of bathymetry change and sediment properties in 

A ^ . They were excluded in the previous model but are included 
in the latest version. 
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b. Case of m=n 



In order to derive an expression for the vector 
mode coupling coefficient for m=n, we differentiate the 
orthonormal condition Eq. (4) using Leibniz rule. The result 
is 



A A 

A^ir.Q) = p^r + Ynn 0 = 



(33) 



1 1 

Pi P 2 



V b H(r,d) ( — - — ) Z, 2 



J ln 



z=H{r ,d) 



Note that, this coupling coefficient is zero for a flat 
horizontal bottom. The latest version of the model has 
included this new term. 



2. SCALAR MODE COUPLING COEFFICIENT, B mn 

Taking the horizontal gradient of both sides of Eq. 
(21), i.e., definition of the vector mode coupling 
coefficient, and applying the Leibniz rule for 
differentiation of a definite integral, we get after some 
manipulations, the following expression for the scalar mode 
coupling coefficient: 
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( 34 ) 



B mn ( r ' 0 ) = \VtA m {r,Q) - fj V h z m .v h z n dz - 



0 



Hi H 2 



z=H(r , 0) 



There is a unique property associated with a complete set of 
orthonormal functions, called the "closure relationship. " For 
the normal modes, this relationship can be expressed as 

T—Z n (z;r,Q) Z n (z';r,Q) =b(z-z') ( 35) 

n p n n 



Taking the horizontal gradient of both sides of the closure 
relationship, multiplying by Z m (z;r,0) and then integrating 
over the entire depth, we get, after some rearranging of 
terms , 



V h Z m {z;r,$) = 



" £ E nn U,Q)Z n (z i r, 6) 



( 36 ) 



where 
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( 37 ) 



Ennir.S) = |^(r,e) + 

W h H(r,d) ( — - — ) Z ln {H-, r ,%) r ,$) 

Pi P 2 



Substituting now Eq. (36) in Eq. (34), we finally obtain the 
following alternative expression for B m : 



= |v^„(r,e) - E E„,(r,e) •E.jlr.e) 



(— — ) Z lm (z; r,6) V h z ln ( z ; r , 6) . V ft H(r, 6) 

Pi P2 



z=H(r,6) 



( 38 ) 



Eq. (38) is valid for both the m?hi and m=n cases. The last 
term of Eq. (38), is new in the model. The magnitude of this 
mode coupling coefficient is generally much less than the 
magnitude of the vector coefficient. 

C. NUMERICAL IMPLEMENTATION 

1 

The major part of the model upgrade was the replacement of 
the old routines with new ones for the computations of the 
rederived mode coupling coefficients according to Eqs . (32), 

(33) and (38) . 

These new routines are contained in a program called 
"sedbot" and are listed in Appendix C. Normal modes and the 
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